A combined density functional theory and molecular mechanics (QM/MM) study of FeCO vibrations in carbonmonoxy myoglobin

Molecular oscillations of ligand motions at the active site of carbonmonoxy myoglobin have been calculated in a protein environment using the combined QM/MM approach. In these calculations, the active site was calculated at a quantum mechanical (QM) level of theory using the B3LYP/6-31 + G* method, while the remaining protein was calculated at the molecular mechanical (MM) level utilizing the AMBER force field. The presence of a torsion mode and the second component of the bending vibration are proof of a bent CO geometry in the heme pocket. 2005 Elsevier B.V. All rights reserved. The heme pocket of myoglobin controls the reactivity of the active site toward O2 and CO adducts [1,2]. Collman et al. [3] proposed that the mechanism for the large difference in affinity between O2 and CO is based on a steric interaction between the ligand adducted to the heme– Fe complex and a distal histidine in close contact with the ligand. In the protein-free heme–Fe–CO complex, a linear orientation of the Fe–C–O bond is preferred, while in myoglobin the CO ligand is bound to the heme–Fe complex in a bent conformation. There are several highresolution X-ray crystallographic structures showing a range of displacements and orientations for the CO ligand in various crystal forms [4–12], and all of them show the heme–Fe–CO geometry to be non-linear. To date, spectroscopic data fail to provide any evidence in support of a non-linear Fe–C–O bond. In the middle-frequency spectral region there is a Fe–CO stretching mode, involving the Fe–C bond, which is experimentally observed [13] at 512 cm . However, there is dispute over the assignment of the Fe–C–O bending modes [14], which are expected to be non-degenerate in the protein due to steric hindrance with the heme pocket. One component of the bending 0009-2614/$ see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2005.11.111 * Corresponding author. Fax: +1 716 645 6505. E-mail address: [email protected] (M. Freindorf). mode is assumed to be at 577 cm 1 in the protein [15]. Unfortunately, there is no experimental confirmation of the second Fe–C–O bending mode, which in the gas phase is almost degenerate with the first component [16]. While there are also expected to be low-frequency Fe–CO bending and torsion modes, their IR and RR intensities are too weak to be observed. Assignment of these low-frequency modes is important because it provides insight into the geometry of the CO ligand bonded to the heme–Fe molecule in the protein. Based on the work presented here, we suggest that previously unassigned weak spectral features in the protein correspond to bending and torsion motions, thereby providing further proof of a bent CO geometry in the heme pocket. In order to assign the low frequency spectral features in carbonmonoxy myoglobin to specific vibrational modes, we carried out a series of combined QM/MM calculations using a B3LYP/6-31 + G*/AMBER potential [17]. The ligand oscillations of the active site in carbonmonoxy myoglobin have already been investigated in the literature by Rovira et al. [18], using a similar QM/MM approach, however the authors report only differences in ligand frequencies between the isolated system in the gas phase and a small QM model system of the protein active site. We report here results of more ambitious calculations which 564 M. Freindorf et al. / Chemical Physics Letters 419 (2006) 563–566 include additional vibrations along with the absolute value of the computed frequencies. Details of the calculations are as follows. A near atomic resolution X-ray structure of Schlichting et al.[11], has the highest resolution (1.15 Å) published to date and accordingly was chosen as the starting point for our study. The 1A6G protein data bank contains three experimental conformers (denoted as A, B, and C) which differ primarily in the position of the distal histidine HIS64. In conformer B, the HIS64 residue is slightly closer to CO than in the A conformer, while in conformer C, the HIS64 residue is swungout of the heme pocket. It is thought that conformers A and B correspond to different protonation states of HIS64, however there is still uncertainty in this assignment. The C conformer has been identified by a diffraction experiment at low pH, which indicates a double protonation state Nde of HIS64. Given the uncertainty associated with the protonation state of HIS64 in the conformers, we calculated both protonation states Nd and Ne of HIS64 in conformers A and B, named as: A-HIS64NdH, A-HIS64NeH, B-HIS64NdH, B-HIS64NeH. Conformer C has been represented by one protonation tautomer C-HIS64NdeH. Hydrogen atoms, which are missing in the experimental structure, were placed in their appropriate positions using the AMBER8 program [19], and their positions were minimized at the molecular mechanical level, utilizing standard parameters for the heme–Fe–CO molecular complex [20] as implemented in the AMBER-8 program. After the energy minimization, each structure was divided into a QM region, which includes the heme–Fe–CO molecule and the proximal histidine HIS93, and a MM region, which includes the remaining protein. The chemical bond between the proximal histidine HIS93 and a side chain of the protein was cut, and resulting free chemical valences were filled out by hydrogen atoms, according to the linking atom method [21,22]. We do not believe the presence of this covalent bond strongly perturbs the CO ligand oscillations because it is located more than 8A from the CO ligand, and the CO ligand oscillations are mostly represented by local Fe–C–O modes. The QM region was simplified by eliminating side chains of the heme ring and a side chain of the HIS93 residue, as was done in previous studies of this system [16,23]. From this starting point, a geometry optimization of the QM system was carried out with the MM system held fixed. The QM/MM calculations included external MM charges in the self consistent field (SCF) computations of the QM wave function, as well as a van der Waals interaction between the QM atoms and the MM atoms. The van der Waals interaction was approximated by a standard Lennard–Jones potential. The QM part of the combined QM/MM calculations was carried out at the B3LYP/6-31 + G* level of theory, which has proven [24] to be very effective for computations of biological systems. The perturbing MM atoms were described by the AMBER-8 force field. The combined QM/MM calculations presented in our study were performed using the Q-Chem [25] program. The geometrical parameters of the heme– Fe–CO active site in the protein and in the gas phase, which were determined in this study, are presented in Table 1. The torsion angle, denoted as a, is defined as the N–Fe– C–O dihedral angle, where the N atom belongs to the heme ring. The bending angle, denoted as b, is the Fe–C–O angle, and the tilting angle, denoted as c, is the N–Fe–C angle, where the N atom belongs to the proximal HIS93 ring. According to our calculations, the optimal value for the torsion angle a is the most sensitive geometric parameter with respect to the various protein conformers; varying from 0 to 100 . However its average value (for the conformers considered in this study) is 31 and is very close to experiment (36 ). The optimal values of the other angles are close to the results of calculations reported by others [23,26] as well as experimental measurements, as indicated in Table 1. The bending angle (b) which was calculated in our study (168 ) is slightly smaller from the value of this angle calculated by Rovira [27] (177 ) and Ryde et al. [23] (171 ), using the similar QM/MM approach. The results of our present calculations are obtained using new values of Lennard–Jones parameters describing the interaction between QM and MM atoms. The new parameters used in our calculations were evaluated using the QM/ MM approach, and they are slightly bigger (about 10– 20%) than parameters used in combined QM/MM calculations by others, which were based on molecular mechanical force fields. Bigger Lennard–Jones parameters cause a stronger repulsive interaction between QM and MM atoms, which causes that the bending angle obtained in our study is slightly smaller than the angle calculated previously by others. We have also performed test calculations for the QM system in the protein having the CO ligand coordinated to the heme iron atom in a linear geometry. The energy difference between the QM system having the optimal bent CO geometry in the protein, and the QM system having the CO ligand bound to the heme molecule in a linear conformation in the protein, was found in our test calculations to have a value between 30 and 50 kcal/mol, depending on the protonation state of distal histidine and depending on the experimental protein conformer (A, B or C). This big CO bending energy which we have estimated in our test calculations in the protein, supports the initial explanation of protein discrimination (lowering affinity of CO binding to heme proteins) suggested many years ago by Collman et al.[3]. Therefore, this issue is a subject of our current investigations. Once the geometry of the QM moiety was optimized, combined QM/MM calculations of the oscillation Hessian were performed, again with the MM portion of the protein held fixed. For comparison, we have also calculated the oscillation Hessian of the QM system in the gas phase (without the protein). The results of our calculations are shown in Table 2, and sketches of vibrational eigenvectors involving large amplitudes of the CO movement are shown in Fig. 1. Three types of vibrational modes were investigated in the present study, namely (1) a low-frequency torsion Table 1 The geometric parameters of the heme–Fe–CO molecule calculated in the protein and in the gas phase, which are compared with the X-ray experimental data Geometrical parameter Fe–C (r) distance (Å) C–O (d) distance (Å) N–Fe–C–O (a) torsion ( ) Fe–C–O (b) bending ( ) N–Fe–C (c) tilting ( ) Gas-phase 1.78 1.15 – 180 180